# -*- coding: utf-8 -*-
import numpy as np
from scipy.constants import e, pi,c, Boltzmann as kb, m_e
from scipy import integrate
import peakdetect as pk
import matplotlib.pyplot as plt
import Graphite_Functions as gf
import FittingFunctions as ff
import time

kb=kb/e  #Boltzmann constant in eV/K
hbar=6.582119624E-16 #Reduced Planck constant in eV*s
az= 3.35E-10 #lattice parameter in z direction. units: m
'''Loading data'''
A1=np.load('Bands data/areaPlus_4.npy') *(2*pi*hbar/az)**2 # Area[E,kz] rows correspond to energy, columns to kz
A2=np.load('Bands data/areaMinus_4.npy') *(2*pi*hbar/az)**2 # Units: eV^2 s^2/m^2 
kz=np.load('Bands data/kz_4.npy') 
Econtours1=np.load('Bands data/Econtours1_4.npy')  
Econtours2=np.load('Bands data/Econtours2_4.npy')
#Fourier transforms
# Pxn1[E,kz,n] with len(n)=120. In the case of
# disconnected energy contours, the first 60 entries in n
# correspond to central energy contour (central pocket)
# and the rest 60 entries to the leg contour.
# In the case of a connected energy contour, the last
# 60 entries are 0
Pxn1=np.load('Bands data/Pxn_plus_6.npy')
Pyn1=np.load('Bands data/Pyn_plus_6.npy') 
w1=2*pi*np.load('Bands data/w_plus_6.npy') #multiplying by 2*pi to get the output of fftfreq in angular frequency. See test2.py
Pxn2=np.load('Bands data/Pxn_minus_6.npy')
Pyn2=np.load('Bands data/Pyn_minus_6.npy')
w2=2*pi*np.load('Bands data/w_minus_6.npy')

'''------------------------- INPUT PARAMETERS -----------------------------'''
#General parameters
Ef_start=-0.024           #Fermi energy in eV
Ef_step=-0.001       #Size of each step between each simulation (Ef_start is being decreased to lower values by this Ef_step value)
ef_count= 1         #number of runs of the simulation, given the Fermi energy decreased by step 

#Main functions
run_sim =  False #Allows the simulation
run_fermi = False #Allows plotting of ds_xx peaks per Ef in dependance on B 
run_plotting = True #Allows plotting of previously acquired data
run_overlap = True #Overlaps measurement data with simulation
run_angle = True #Allows plotting of angle dependance on B for measurement data
run_angle_comparison = True #Allows plotting of angle dependance on B for all measurement data in one graph

#Fermi energy plotting custom parameters
delete_peaks_fermi = False #deletes peaks that are not useful for fermi energy calculation (only 3 peaks needed)
deleted_peaks = 1 #number of peaks to delete

#Plotting, if run_plotting is true
plot_sxx_to_B = False #Plots sigma_xx in dependance on B
plot_dsxx_to_B = True #Plots derigation of sigma_xx in dependance on B
plot_sxx_to_wwc = False #plots sigma_xx in dependance on omega/omega c0
plot_dsxx_to_wwc = False #plots derivation of sigma_xx in dependance on omega/omega c0

#Plotting, plot_angle_sum must be true to enable plot_Bperpendicular and plot_volt_dependance
#works properly for Fr_det_min = 0 and Fr_det_max = 0.425
plot_angle_onebyone = False #Plots every single measurement data from grenoble, one by one, with Fr_det_min+max resolution
plot_angle_sum = True #Plots only the peaks development with angle in dependance on Bexp
plot_Bperpendicular = True #Plots only the peaks development with angle in dependance on perpendicular component of Bexp
#plot_volt_dependance is for config only, often causes error because of different length of lists of detected peaks, which is dependant on the selected Fr_sum_min/max field window because of normalization of these data
plot_volt_dependance = True #Plots voltage dependance of data to ensure that zero angle value is correct
custom_peak_detection = True #Detect missing peaks for Fr_det_min = 0 and Fr_det_max = 0.425 in grenoble file 2 and 6
#might be necessary to set custom_peak_delete to True if running plot_Bperpendicular and plot_volt_dependance
custom_peak_delete = True #Deletes peaks for Fr_det_min = 0 and Fr_det_max = 0.425 that emerge with higher magnetic field from file 13 up
smooth_const = 25 #25, smoothing constant of grenoble measurement data
lookahead_const = 15 #15, lookahead constant of grenoble measurement data for peak detection
peakdetect_const = 0.15 #0.15, minimum value to detect peaks in grenoble measurement data
counter_angle_original = 2 #2, first measurement data file from grenoble with usable peaks
last_measurement_data = 18 #18, last measurement data file from grenoble with usable peaks
multiplicator_angle_comparison = 2.85 #multiplicates the resulting signal so its easier to see the peaks

#Overlap of measurement data with simulation data
plot_complete_view = True #Plots complete view of overlap of measurement data with simulation data
plot_detailed_view = True #Plots detailed view of overlap of measurement data with simulation data
Fr_det_min = 0.01 #detailed section of that field, min value
Fr_det_max = 0.425 #detailed section of that field, max value
Fr_sum_min = 0 #detailed section of that field, min value for summary of run_angle
Fr_sum_max = 0.425 #detailed section of that field, max value for summary of angles run_angle
Fr_comp_min = -0.425 #detailed section of field for run_angle_comparison, min value
Fr_comp_max = 0.425 #detailed section of field for run_angle_comparison, max value
skiprowsvar = 1 #30 for classical files, 1 for corrected files
MeasDataFileName="Experimental data/Bulgarian graphite HFEPR/corrected_1-Graphine Bulgarian_2--2T 210GHz.dat" #Path to measurement data
data_offset = 0.0266 #0.0247 #zero is shifted by 24,7mT
Fr_val = 16.5 #Magnetic field range in the plot, must be high enough
meas_full_y_offset = 0 #0 for normal data, -0.905 for zero corrected data
meas_full_multiplier = 1 #1 for normal data, 10 for zero corrected data
meas_det_y_offset = 0#0 for normal data, -0.94 for zero corrected data
meas_det_multiplier = 1 #1 for normal data, 14 for zero corrected data

#Simulation parameters
B=np.append(np.arange(0.01,0.15,0.0002), np.arange(0.15,0.425,0.001))
f=190E9#measurement frequency of microwaves
T=5 #temperature in K
a=0.08E-3 #proportionality constant when G depends on sqrt(B):  G=a*sqrt(B). a in eV
G= 20E-6 #broadening parameter in eV, now completelyreplaced with a
mc0=0.06*m_e #effective mass of the carriers in graphite, in kg

#Peak detection parameters, normalized
delta_s_b = 0.01 #min. difference in sigma in dependance on B to be recognized as peak
delta_ds_b = 0.01 #min. diff in deriv. of sigma in dependance on B
delta_s_wwc = 0.01 #min. difference in sigma in dependance on w/wc0 to be recognized as peak
delta_ds_wwc = 0.01 #min. difference in derivation of sigma in dependance on w/wc0 to be recognized as peak
la_range = 30 #Lookahead range, how many points is compared between each other

w=2*pi*f #Microwave frequency in Hz
wc0= e*B/mc0 #cyclotron frequency w.r.t the magnetic field. in Hz
ef_counter = 0
Ef=Ef_start

'''Calculating the optical conductivity'''
def conductivity3(B,w,a):
    '''Calculates the optical conductivity, accepts B as a vector. Now w1 and w2 are real cyclotron frequencies 
    instead of integer numbers'''
    t1=time.time()
    sigma_xx=np.zeros(len(B))
    dE1=Econtours1[1]-Econtours1[0]
    dE2=Econtours2[1]-Econtours2[0]
    dkz=kz[1]-kz[0]
    for ib,BB in enumerate(B):
        print('magnetic field', ib, ' out of ', len(B) )
        dummy1_cent= np.nansum(n1_cent**2 *(Py1_cent/BB)**2/( hbar**2 *(w -BB*w1_cent )**2 +(a*np.sqrt(BB))**2), -1)
        dummy1_leg= np.nansum(n1_leg**2 *(Py1_leg/BB)**2/( hbar**2 *(w -BB*w1_leg )**2 +(a*np.sqrt(BB))**2), -1) 
        dummy2_cent= np.nansum(n2_cent**2 *(Py2_cent/BB)**2/( hbar**2 *(w -BB*w2_cent )**2 +(a*np.sqrt(BB))**2), -1) 
        dummy2_leg= np.nansum(n2_leg**2 *(Py2_leg/BB)**2/( hbar**2 *(w -BB*w2_leg )**2 +(a*np.sqrt(BB))**2), -1)  
        sigma_xx[ib]= np.nansum(-1*df1[:,np.newaxis]*(a*np.sqrt(BB))*(m1_cent/m1_cent**2) *dummy1_cent)*dE1*dkz*2*pi/az \
        +np.nansum(-1*df1[:,np.newaxis]*3*(a*np.sqrt(BB))*(m1_leg/m1_leg**2) *dummy1_leg)*dE1*dkz*2*pi/az \
        +np.nansum(-1*df2[:,np.newaxis]*(a*np.sqrt(BB))*(m2_cent/m2_cent**2) *dummy2_cent)*dE2*dkz*2*pi/az  \
        +np.nansum(-1*df2[:,np.newaxis]*3*(a*np.sqrt(BB))*(m2_leg/m2_leg**2) *dummy2_leg)*dE2*dkz*2*pi/az
        print('sigma_xx[ib]', sigma_xx[ib])
    t2=time.time()
    print('time for conductivity calculation with for loop ', t2-t1)
    return sigma_xx

'''MAIN LOOP'''
if run_sim == True:
    while ef_counter < ef_count:  
        '''calculating df/dE to simplify calculations later'''
        df1= gf.dfdE(Econtours1, Ef, T)
        df2= gf.dfdE(Econtours2, Ef, T)
        
        '''slicing Econtour so we perform less calculations'''
        print('length before slicing E1, E2', len(Econtours1), len(Econtours2))
        mask1=np.where(1E-1<np.abs(df1))
        mask2=np.where(1E-1<np.abs(df2))
        
        Econtours1=Econtours1[mask1]
        Econtours2=Econtours2[mask2]
        A1=A1[mask1[0],:]
        A2=A2[mask2[0],:]
        Pxn1=Pxn1[mask1[0],:,:]
        Pyn1=Pyn1[mask1[0],:,:]
        w1=w1[mask1[0],:,:]
        Pxn2=Pxn2[mask2[0],:,:]
        Pyn2=Pyn2[mask2[0],:,:]
        w2=w2[mask2[0],:,:]
        
        df1= gf.dfdE(Econtours1, Ef, T)
        df2= gf.dfdE(Econtours2, Ef, T)
        print('length after slicing E1, E2', len(Econtours1), len(Econtours2))
            
        '''Separating components of cyclotron orbits'''
        Px1_cent, Px1_leg= Pxn1[:,:,:60], Pxn1[:,:,60:]
        Py1_cent, Py1_leg= Pyn1[:,:,:60], Pyn1[:,:,60:]
        w1_cent, w1_leg= w1[:,:,:60], w1[:,:,60:]
        Px2_cent, Px2_leg= Pxn2[:,:,:60], Pxn2[:,:,60:]
        Py2_cent, Py2_leg= Pyn2[:,:,:60], Pyn2[:,:,60:]
        w2_cent, w2_leg= w2[:,:,:60], w2[:,:,60:]
        
        '''dealing with problematic points: they are close to the points where both bands
        touch each other. The Fourier transform is inaccurate in such points, that is why
        we turn them into zero'''
        #The bands touch at four points, the center point and three legs
        ind1_cent= np.where(Py1_cent[:,:,1]!=np.max(Py1_cent,-1))
        ind1_leg= np.where(Py1_leg[:,:,1]!=np.max(Py1_leg,-1))
        ind2_cent= np.where(Py2_cent[:,:,1]!=np.max(Py2_cent,-1))
        ind2_leg= np.where(Py2_leg[:,:,1]!=np.max(Py2_leg,-1))
        Px1_cent[ind1_cent[0],ind1_cent[1],:] = np.zeros(60) 
        Py1_cent[ind1_cent[0],ind1_cent[1],:] = np.zeros(60)
        Px1_leg[ind1_leg[0],ind1_leg[1],:] = np.zeros(60)
        Py1_leg[ind1_leg[0],ind1_leg[1],:] = np.zeros(60)
        Px2_cent[ind2_cent[0],ind2_cent[1],:] = np.zeros(60)
        Py2_cent[ind2_cent[0],ind2_cent[1],:] = np.zeros(60)
        Px2_leg[ind2_leg[0],ind2_leg[1],:] = np.zeros(60)
        Py2_leg[ind2_leg[0],ind2_leg[1],:] = np.zeros(60)
        
        '''frequencies translated to integer numbers: easier for optical 
        conductivity calculations'''
        n1_cent, n1_leg= w1_cent/w1_cent[:,:,1,np.newaxis], w1_leg/w1_leg[:,:,1,np.newaxis]
        n2_cent, n2_leg= w2_cent/w2_cent[:,:,1,np.newaxis], w2_leg/w2_leg[:,:,1,np.newaxis]
        n1_cent[np.isnan(n1_cent)], n1_leg[np.isnan(n1_leg)] = 0, 0
        n2_cent[np.isnan(n2_cent)], n2_leg[np.isnan(n2_leg)] = 0, 0
        
        '''Calculating the effective mass'''
        # Effective mass in eV*s^2/m^2. Calculated as eB/wc (in S.I. units,
        # divided by e to get it in eV) B=1T in Fourier transf. calculation 
        m1_cent, m1_leg =e*1/w1_cent[:,:,1]/e, e*1/w1_leg[:,:,1]/e 
        m2_cent, m2_leg =-e*1/w2_cent[:,:,1]/e, -e*1/w2_leg[:,:,1]/e 
        m1_cent[np.isinf(m1_cent)], m1_leg[np.isinf(m1_leg)] = 0, 0
        m2_cent[np.isinf(m2_cent)], m2_leg[np.isinf(m2_leg)] = 0, 0
        print('m2_center', m2_cent)
        
        sigma_xx=conductivity3(B,w,a)
        ds_xx=np.gradient(sigma_xx, B[1]-B[0])
        
        '''Saving data'''
        file_title = 'Ef_Sim/Ef_' +  '{:.1f}'.format(Ef*1000) + 'meV.npz'
        np.savez(file_title, B=B, sigma_xx=sigma_xx, ds_xx=ds_xx, wwc0=w/wc0)
    
        Ef = Ef + Ef_step
        ef_counter += 1

'''Plotting of data, displays all the data of every simulation'''
if run_plotting == True:
    ef_counter = 0
    Ef = Ef_start
    while ef_counter < ef_count:
        name= '$f$ = '+'{:n}'.format(f*1E-9) + ' GHz, $T$ = ' +'{:n}'.format(T) + r' K, $E_{F}$ = '+'{:n}'.format(Ef*1000)+' meV, $a$ = '+str(a*1000)+' meV'
        file_title = 'Ef_Sim/Ef_' +  '{:.1f}'.format(Ef*1000) + 'meV.npz'
        Loaded_data = np.load(file_title)   
        '''Data normalization'''
        sigma_normalized = []
        ds_xx_normalized = []
        sval_max = abs(max(Loaded_data['sigma_xx'], key=abs))
        dsval_max = abs(max(Loaded_data['ds_xx'], key=abs))
        for svar in Loaded_data['sigma_xx']:
            sigma_normalized.append(svar/sval_max)
        for dsvar in Loaded_data['ds_xx']:
            ds_xx_normalized.append(dsvar/dsval_max)
        ''' PLOTS'''
        if  plot_sxx_to_B == True:        
            fig1_peaks = pk.peakdetect(sigma_normalized, Loaded_data['B'], lookahead=la_range, delta=delta_s_b) 
            fig1_hPeaks = np.array(fig1_peaks[0])
            fig1_lPeaks = np.array(fig1_peaks[1])
            fig,ax=plt.subplots()
            ax.plot(Loaded_data['B'], sigma_normalized)
            ax.plot(fig1_hPeaks[:,0], fig1_hPeaks[:,1], 'r+')
            plt.axis([Fr_det_min, Fr_det_max, -0.1, 1.1])
            plt.xlabel('$B$ (T)', fontsize=15)            
            plt.ylabel('$\sigma_{xx}$ (arb.u.)', fontsize=15)
            plt.title(name)                          
        if plot_dsxx_to_B == True:
            fig2_peaks = pk.peakdetect(ds_xx_normalized, Loaded_data['B'], lookahead=la_range, delta=delta_ds_b)
            fig2_hPeaks = np.array(fig2_peaks[0])
            fig2_lPeaks = np.array(fig2_peaks[1])
            fig,ax=plt.subplots()
            ax.plot(Loaded_data['B'], ds_xx_normalized)
            ax.plot(fig2_hPeaks[:,0], fig2_hPeaks[:,1], 'r+')
            ax.plot(fig2_lPeaks[:,0], fig2_lPeaks[:,1], 'k+')
            plt.axis([Fr_det_min, Fr_det_max, -1.1, 1.1])
            plt.xlabel('$B$ (T)', fontsize=18)          
            plt.xticks(fontsize=13)
            plt.yticks(fontsize=13)
            plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=18)  
            plt.title(name, fontsize=17)  
            plt.gcf().subplots_adjust(left=0.16)
            plt.gcf().subplots_adjust(bottom=0.14)
        if  plot_sxx_to_wwc == True:
            fig3_peaks = pk.peakdetect(sigma_normalized, Loaded_data['wwc0'], lookahead=la_range, delta=delta_s_wwc) 
            fig3_hPeaks = np.array(fig3_peaks[0])
            fig3_lPeaks = np.array(fig3_peaks[1])
            fig,ax=plt.subplots()
            ax.plot(Loaded_data['wwc0'],sigma_normalized)
            plt.axis([0,20, -0.1, 1.1])
            plt.xticks([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])
            plt.xlabel('$\omega/\omega_c$', fontsize=15)            
            plt.ylabel('$\sigma_{xx}$ (arb.u.)', fontsize=15) 
            plt.title(name)  
        if plot_dsxx_to_wwc == True:
            fig4_peaks = pk.peakdetect(ds_xx_normalized, Loaded_data['wwc0'], lookahead=la_range, delta=delta_ds_wwc) 
            fig4_hPeaks = np.array(fig4_peaks[0])
            fig4_lPeaks = np.array(fig4_peaks[1])
            fig,ax=plt.subplots()
            ax.plot(Loaded_data['wwc0'],ds_xx_normalized)
            plt.axis([0,20, -1.1, 1.1])
            plt.xlabel('$\omega/\omega_c$', fontsize=15)            
            plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=15) 
            plt.title(name)  
        del Loaded_data.f
        Loaded_data.close()
        Ef = Ef + Ef_step
        ef_counter += 1

'''Plotting of data, displays only the dependance of peaks on Fermi energy and magnetic field'''
if run_fermi == True:
    ef_counter = 0
    Ef = Ef_start
    datalist = []
    sizelist = []
    pk_ef = np.array([])
    fig,ax=plt.subplots()
    plt.axis([Fr_det_min,Fr_det_max, Ef_start*1000+(Ef_step*1000*ef_count)-(Ef_step*1000/2), Ef_start*1000-(Ef_step*1000/2)])
    plt.xlabel('$B$ (T)', fontsize=15)            
    plt.ylabel('$E_{F}$ (meV)', fontsize=15) 
    while ef_counter < ef_count:
        file_title = 'Ef_Sim/Ef_' +  '{:.1f}'.format(Ef*1000) + 'meV.npz'
        Loaded_data = np.load(file_title)
        
        '''Data normalization'''
        ds_xx_normalized = []
        dsval_max = abs(max(Loaded_data['ds_xx'], key=abs))
        for dsvar in Loaded_data['ds_xx']:
            ds_xx_normalized.append(dsvar/dsval_max)  
        s_xx_normalized = []
        sval_max = abs(max(Loaded_data['sigma_xx'], key=abs))
        for svar in Loaded_data['sigma_xx']:
            s_xx_normalized.append(svar/sval_max)    
        ef_peaks = pk.peakdetect(s_xx_normalized, Loaded_data['B'], lookahead=la_range, delta=delta_ds_b)
        ef_hPeaks = np.array(ef_peaks[0])
        if delete_peaks_fermi == True:
            if Ef < -0.02:
                ef_hPeaks = np.delete(ef_hPeaks, 0, 0)
            counter_aux = 0
            while counter_aux <= (deleted_peaks-1):
                ef_hPeaks = np.delete(ef_hPeaks, 0, 0)
                counter_aux += 1
            ef_hPeaks = np.delete(ef_hPeaks, 4, 0)
            
        ax.plot(ef_hPeaks[:,0], np.full(len(ef_hPeaks[:,1]), Ef*1000), 'ro')
        print("Shape of peaks for Ef= "+'{:.1f}'.format((Ef)*1000)+": ")
        print(ef_hPeaks.shape)
        datalist.append(ef_hPeaks[:,0])
        sizelist.append(ef_hPeaks[:,0].shape)
        
        del Loaded_data.f
        Loaded_data.close()
        Ef = Ef + Ef_step
        ef_counter += 1
    plt.title('$f$ = '+'{:n}'.format(f*1E-9) + ' GHz, $T$ = ' +'{:n}'.format(T) + r' K, $E_{F}$ from '+'{:.1f}'.format(Ef_start*1000)+' to '+ '{:.1f}'.format((Ef-Ef_step)*1000)+' meV')    
    
    print("List of all detected peaks")
    print(datalist)
    print("List of shapes of these peaks")
    print(sizelist)
    print("Largest value on previous list")
    largest_val=max(sizelist)
    print(largest_val[0])
    Ef = Ef_start 
    
    pk_1 = np.array([])
    pk_2 = np.array([])
    pk_3 = np.array([])  
    pk_4 = np.array([]) 
    for sup_x in range(ef_count):
        for sup_y in range(largest_val[0]):            
            if sup_y == 0:
                pk_1 = np.append(pk_1,datalist[sup_x][sup_y])
            if sup_y == 1:
                pk_2 = np.append(pk_2,datalist[sup_x][sup_y])
            if sup_y == 2:
                pk_3 = np.append(pk_3,datalist[sup_x][sup_y])   
            if sup_y == 3:
                pk_4 = np.append(pk_4,datalist[sup_x][sup_y])   
        pk_ef = np.append(pk_ef,Ef*1000)        
        Ef = Ef + Ef_step
        
    print("Peaks sorted by fermi energy")       
    print(pk_1)
    print(pk_2)
    print(pk_3)
    print(pk_4)
    print("Corresponding Fermi energy list of the same size")
    print(pk_ef)
    
    eq1 =  np.polyfit(pk_1, pk_ef, 2)
    eq2 =  np.polyfit(pk_2, pk_ef, 2)
    eq3 =  np.polyfit(pk_3, pk_ef, 2)
    eq4 =  np.polyfit(pk_4, pk_ef, 2)
    ax.plot(pk_1, eq1[0]*pk_1**2+eq1[1]*pk_1+eq1[2],'g')
    ax.plot(pk_2, eq2[0]*pk_2**2+eq2[1]*pk_2+eq2[2],'b')
    ax.plot(pk_3, eq3[0]*pk_3**2+eq3[1]*pk_3+eq3[2],'y')
    ax.plot(pk_4, eq4[0]*pk_4**2+eq4[1]*pk_4+eq4[2],'k')
            
    print("Equation of the first peaks set (green): " +str(eq1[0])+" * x^2 + "+str(eq1[1])+" * x + "+str(eq1[2]))
    print("Equation of the second peaks set (blue): " +str(eq2[0])+" * x^2 + "+str(eq2[1])+" * x + "+str(eq2[2]))
    print("Equation of the third peaks set (yellow): " +str(eq3[0])+" * x^2 + "+str(eq3[1])+" * x + "+str(eq3[2]))
    print("Equation of the fourth peaks set (black): " +str(eq4[0])+" * x^2 + "+str(eq4[1])+" * x + "+str(eq4[2]))

'''Measurement data overlaping with simulation'''
if run_overlap == True:
    ef_counter = 0
    Ef = Ef_start
 
    MeasData = np.loadtxt(open(MeasDataFileName,'r'),skiprows=skiprowsvar)
    B_exp=(MeasData[:,0]/1000)+data_offset
    ds_exp=MeasData[:,1]
    B_exp_det = B_exp
    ds_exp_det = ds_exp
    B_exp = B_exp[np.logical_and(-Fr_val<B_exp,B_exp<Fr_val)]
    B_exp_det = B_exp_det[np.logical_and(Fr_det_min<B_exp_det,B_exp_det<Fr_det_max)]
    ds_exp = 1E6*ds_exp[np.logical_and(-Fr_val<B_exp,B_exp<Fr_val)]
    ds_exp_det = ds_exp[np.logical_and(Fr_det_min<B_exp,B_exp<Fr_det_max)]
    
    '''Measurement data normalization'''
    meas_full_normalized = []
    meas_detail_normalized = []
    meas_full_max = abs(max(ds_exp, key=abs))
    meas_detail_max = abs(max(ds_exp_det, key=abs))
    for mvar in ds_exp:
        meas_full_normalized.append(((mvar/meas_full_max)*meas_full_multiplier)+meas_full_y_offset)
    for nvar in ds_exp_det:
        meas_detail_normalized.append(((nvar/meas_detail_max)*meas_det_multiplier)+meas_det_y_offset)
    
    while ef_counter < ef_count:
        name= '{:n}'.format(f*1E-9) + 'GHz_'+'{:n}'.format(Ef*1000) + 'meV_' + '{:n}'.format(T) + 'K'
        file_title = 'Ef_Sim/Ef_' +  '{:.1f}'.format(Ef*1000) + 'meV.npz'
        Loaded_data = np.load(file_title)
        
        '''Simulation data normalization'''
        ds_xx_normalized = []
        dsval_max = abs(max(Loaded_data['ds_xx'], key=abs))
        for dsvar in Loaded_data['ds_xx']:
            ds_xx_normalized.append(dsvar/dsval_max)
        s_xx_normalized = []
        sval_max = abs(max(Loaded_data['sigma_xx'], key=abs))
        for svar in Loaded_data['sigma_xx']:
            s_xx_normalized.append(svar/sval_max)
        
        if plot_complete_view == True:
            fig,ax=plt.subplots() #overlap of measurement with simulation, complete view
            ax.plot(B_exp, meas_full_normalized, 'b')
            plt.xlabel('$B_{exp} (T)$', fontsize=15)            
            plt.ylabel('$d\sigma_{xx}/dB$', fontsize=15)  
            plt.title("$E_{f}$ = "+ '{:.1f}'.format((Ef)*1000)+'meV')
        
        if plot_detailed_view == True:
            fig,ax=plt.subplots() #overlap of measurement with simulation, detailed view
            ax.plot(B_exp_det, meas_detail_normalized, 'b')
            ax.plot(Loaded_data['B'], ds_xx_normalized, 'r')
            plt.axis([Fr_det_min, Fr_det_max, -1.1, 1.1])
            plt.xlabel('$B$ (T)', fontsize=20)            
            plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=20)  
            plt.title("$E_{f}$ = "+ '{:.1f}'.format((Ef)*1000)+' meV, $a$ = '+str(a*1000)+' meV', fontsize=20)
            plt.xticks(fontsize = 13)
            plt.yticks(fontsize = 13)
            plt.gcf().subplots_adjust(left=0.16)
            plt.gcf().subplots_adjust(bottom=0.14)
                  
        integrated_s = integrate.cumtrapz(ds_exp_det,B_exp_det)
        B_exp_det2 = B_exp_det[1:]
        lin_params=np.polyfit(B_exp_det2,integrated_s,1)
        base= lin_params[1] + lin_params[0]*B_exp_det2
        integrated_s = integrated_s - base   
        ints_xx_normalized = []
        intsval_max = abs(max(integrated_s, key=abs))
        for intsvar in integrated_s:
            ints_xx_normalized.append(((intsvar/intsval_max)*0.8)+0.15)
            
        fig,ax=plt.subplots() #overlap of measurement with simulation, detailed view
        ax.plot(B_exp_det2, ints_xx_normalized, 'b')
        ax.plot(Loaded_data['B'], s_xx_normalized, 'r')
        plt.axis([Fr_det_min, Fr_det_max, -1.1, 1.1])
        plt.xlabel('$B (T)$', fontsize=15)            
        plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=15)  
        plt.title("$E_{f}$ = "+ '{:.1f}'.format((Ef)*1000)+'meV')
            
        del Loaded_data.f
        Loaded_data.close()
        Ef = Ef + Ef_step
        ef_counter += 1
    
def load_angle_data(dummy_counter):
    switcher={
        1: "Experimental data/Rotated_graphite/001-190GHz-NG-0-401V.dat",
        2: "Experimental data/Rotated_graphite/002-190GHz-NG-0-453V.dat",
        3: "Experimental data/Rotated_graphite/003-190GHz-NG-0-502V.dat",
        4: "Experimental data/Rotated_graphite/004-190GHz-NG-0-552V.dat",
        5: "Experimental data/Rotated_graphite/005-190GHz-NG-0-602V.dat",
        6: "Experimental data/Rotated_graphite/006-190GHz-NG-0-651V.dat",
        7: "Experimental data/Rotated_graphite/007-190GHz-NG-0-700V.dat",
        8: "Experimental data/Rotated_graphite/008-190GHz-NG-0-752V.dat",
        9: "Experimental data/Rotated_graphite/009-190GHz-NG-0-800V.dat",
        10: "Experimental data/Rotated_graphite/010-190GHz-NG-0-852V.dat",
        11: "Experimental data/Rotated_graphite/011-190GHz-NG-0-902V.dat",
        12: "Experimental data/Rotated_graphite/012-190GHz-NG-0-954V.dat",
        13: "Experimental data/Rotated_graphite/013-190GHz-NG-1-002V.dat",
        14: "Experimental data/Rotated_graphite/014-190GHz-NG-1-050V.dat",
        15: "Experimental data/Rotated_graphite/015-190GHz-NG-1-103V.dat",
        16: "Experimental data/Rotated_graphite/016-190GHz-NG-1-151V.dat",
        17: "Experimental data/Rotated_graphite/017-190GHz-NG-1-205V.dat",
        18: "Experimental data/Rotated_graphite/018-190GHz-NG-1-250V.dat"
        }
    return switcher.get(dummy_counter, "Error: File does not exist")
    
def load_angle(dummy_counter):
    #var_angle = -36.14
    var_angle = -35.4
    switcher={
        1: ((0.401/0.01389)+var_angle),
        2: ((0.453/0.01389)+var_angle),
        3: ((0.502/0.01389)+var_angle), #0°
        4: ((0.552/0.01389)+var_angle),
        5: ((0.602/0.01389)+var_angle), 
        6: ((0.651/0.01389)+var_angle), 
        7: ((0.700/0.01389)+var_angle),
        8: ((0.752/0.01389)+var_angle),
        9: ((0.800/0.01389)+var_angle),
        10: ((0.852/0.01389)+var_angle),
        11: ((0.902/0.01389)+var_angle),
        12: ((0.954/0.01389)+var_angle),
        13: ((1.002/0.01389)+var_angle),
        14: ((1.050/0.01389)+var_angle),
        15: ((1.103/0.01389)+var_angle),
        16: ((1.151/0.01389)+var_angle),
        17: ((1.205/0.01389)+var_angle),
        18: ((1.250/0.01389)+var_angle)
        }
    return switcher.get(dummy_counter, "Error: Angle does not exist")

def load_offset(dummy_counter):
    #'''
    switcher={
        1: 0,
        2: 0.003026,
        3: -0.0055945,
        4: -0.0053485,
        5: 0.004586,
        6: -0.005935,
        7: 0.0047735,
        8: -0.00519,
        9: 0.0052195,
        10: -0.005102,
        11: 0.0052605,
        12: -0.00519,
        13: 0.005657,
        14: -0.0061285,
        15: 0.0071605,
        16: -0.006703,
        17: 0.00563,
        18: -0.0081165
        }
    return switcher.get(dummy_counter, "Error: Offset does not exist")    

'''Angle depentance on magnetic field plotting'''
if run_angle == True:

    if plot_angle_sum == True:
        fig2, ax2=plt.subplots()
        plt.xlabel('$B$ (T)', fontsize=20)            
        plt.ylabel('$Úhel$ (°)', fontsize=20)  
        plt.title("Vývoj vybraných vrcholů", fontsize=20)
        plt.xticks(fontsize=15)
        plt.yticks(fontsize=15)
        plt.gcf().subplots_adjust(bottom=0.15)
        
        fig6, ax6=plt.subplots()
        plt.xlabel('$B$ (T)', fontsize=20)            
        plt.ylabel('$Úhel$ (°)', fontsize=20)  
        plt.title("Vývoj všech vrcholů", fontsize=20)
        plt.xticks(fontsize=15)
        plt.yticks(fontsize=15)
        plt.gcf().subplots_adjust(bottom=0.15)
        
        if  plot_Bperpendicular == True:
            fig3, ax3=plt.subplots()
            plt.xlabel(r'$B_{\bot}$ (T)', fontsize=15)            
            plt.ylabel('$Úhel$ (°)', fontsize=15)  
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
            plt.title("Úhel od "+'{:.1f}'.format(load_angle(counter_angle_original))+"° do "+'{:.1f}'.format(load_angle(last_measurement_data))+"°, kosinus $\Phi$", fontsize=15)
        if plot_volt_dependance == True:    
            fig4, ax4=plt.subplots()
            plt.xlabel('$1/B$ (T$^{-1}$)', fontsize=15)            
            plt.ylabel(r'cos($\theta$)', fontsize=15)  
            plt.title("Vývoj napětí", fontsize=15)
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
        
    datalist_angle = []
    sizelist_angle = []
    list_x = []
    list_y = []
    pk_angle = np.array([])
    counter_angle = counter_angle_original
    while counter_angle <= last_measurement_data:
        loaded_angle_data = load_angle_data(counter_angle)
        AngleData = np.loadtxt(open(loaded_angle_data,'r'),skiprows=30) 
        loaded_offset = load_offset(counter_angle)
        
        B_ang = (AngleData[:,0])+loaded_offset
        ds_ang = AngleData[:,1]
        B_ang_det = B_ang
        ds_ang_det = ds_ang
        
        B_ang = B_ang[np.logical_and(-Fr_val<B_ang,B_ang<Fr_val)]
        B_ang_det = B_ang_det[np.logical_and(Fr_sum_min<B_ang_det,B_ang_det<Fr_sum_max)]
        ds_ang = ds_ang[np.logical_and(-Fr_val<B_ang,B_ang<Fr_val)]
        ds_ang_det = ds_ang[np.logical_and(Fr_sum_min<B_ang,B_ang<Fr_sum_max)]     
        ds_smooth = ff.smooth(ds_ang_det, smooth_const)
    
        '''Measurement data normalization'''
        angle_normalized = []
        angle_max = abs(max(ds_smooth, key=abs))
        for ovar in ds_smooth:
            angle_normalized.append(ovar/angle_max)
         
        '''Peak Detect'''
        angle_peaks = pk.peakdetect(angle_normalized, B_ang_det, lookahead=lookahead_const, delta=peakdetect_const) 
        angle_hPeaks = np.array(angle_peaks[0])
        angle_lPeaks = np.array(angle_peaks[1])
        
        if custom_peak_detection == True:
            if ((counter_angle == 3) and (Fr_det_min <= 0.185) and (Fr_det_max >= 0.2015)) == True:
                B_ang_custom3 = B_ang
                B_ang_custom3 = B_ang_custom3[np.logical_and(0.185<B_ang_custom3,B_ang_custom3<0.2015)] #0.2015 first peak, 0.2125 second peak
                ds_ang_custom3 = ds_ang[np.logical_and(0.185<B_ang,B_ang<0.2015)]
                ds_smooth_custom3 = ff.smooth(ds_ang_custom3, smooth_const)
                angle_normalized_custom3 = []
                angle_max_custom3 = abs(max(ds_smooth, key=abs))
                for ovar_custom3 in ds_smooth_custom3:
                    angle_normalized_custom3.append(ovar_custom3/angle_max_custom3)
                angle_peaks_custom3 = pk.peakdetect(angle_normalized_custom3, B_ang_custom3, lookahead=10, delta=0.001)
                angle_hPeaks_custom3 = np.array(angle_peaks_custom3[0])
                angle_hPeaks = np.concatenate((angle_hPeaks, angle_hPeaks_custom3), axis=0)
                
            if ((counter_angle == 7) and (Fr_det_min <= 0.0275) and (Fr_det_max >= 0.0335)) == True:
                B_ang_custom7 = B_ang
                B_ang_custom7 = B_ang_custom7[np.logical_and(0.0275<B_ang_custom7,B_ang_custom7<0.0335)]
                ds_ang_custom7 = ds_ang[np.logical_and(0.0275<B_ang,B_ang<0.0335)]
                ds_smooth_custom7 = ff.smooth(ds_ang_custom7, smooth_const)
                angle_normalized_custom7 = []
                angle_max_custom7 = abs(max(ds_smooth, key=abs))
                for ovar_custom7 in ds_smooth_custom7:
                    angle_normalized_custom7.append(ovar_custom7/angle_max_custom7)
                angle_peaks_custom7 = pk.peakdetect(angle_normalized_custom7, B_ang_custom7, lookahead=10, delta=0.01)
                angle_hPeaks_custom7 = np.array(angle_peaks_custom7[0])
                angle_hPeaks = np.concatenate((angle_hPeaks, angle_hPeaks_custom7), axis=0)
        
        '''Plotting'''
        if plot_angle_onebyone == True:
            x_hpeaks_plus = []
            x_lpeaks_plus = []
            x_hpeaks_minus = []
            x_lpeaks_minus = []
            fig,ax=plt.subplots() #angle data plotting
            ax.plot(B_ang_det, angle_normalized, 'b')
            largest_val_hpeaks = max(angle_hPeaks.shape)
            largest_val_lpeaks = max(angle_lPeaks.shape)
            for y_hvar in range(largest_val_hpeaks):  
                x_hpeaks_actual = angle_hPeaks[:,0][y_hvar]
                y_hpeaks_actual = angle_hPeaks[:,1][y_hvar]
                if x_hpeaks_actual <= 0:
                    ax.plot(x_hpeaks_actual, y_hpeaks_actual, 'k+')
                    x_hpeaks_minus.append(x_hpeaks_actual)
                if x_hpeaks_actual >= 0:
                    ax.plot(x_hpeaks_actual, y_hpeaks_actual, 'r+')
                    x_hpeaks_plus.append(x_hpeaks_actual)
            for y_lvar in range(largest_val_lpeaks):  
                x_lpeaks_actual = angle_lPeaks[:,0][y_lvar]
                y_lpeaks_actual = angle_lPeaks[:,1][y_lvar]
                if x_lpeaks_actual <= 0:
                    ax.plot(x_lpeaks_actual, y_lpeaks_actual, 'g+')
                    x_lpeaks_minus.append(x_lpeaks_actual)
                if x_lpeaks_actual >= 0:
                    ax.plot(x_lpeaks_actual, y_lpeaks_actual, 'y+')
                    x_lpeaks_plus.append(x_lpeaks_actual)
            x_lpeaks_minus=np.sort(x_lpeaks_minus)
            x_hpeaks_minus=np.sort(x_hpeaks_minus)
            x_lpeaks_plus=np.sort(x_lpeaks_plus)
            x_hpeaks_plus=np.sort(x_hpeaks_plus)
            print(" ")
            print("File MD"+str(counter_angle)+", angle: "+'{:.1f}'.format(load_angle(counter_angle))+"°")
            print("One-by-one detected black peaks:"+str(x_hpeaks_minus))
            print("One-by-one detected green peaks: "+str(x_lpeaks_minus))
            print("One-by-one detected red peaks: "+str(x_hpeaks_plus))
            print("One-by-one detected yellow peaks: "+str(x_lpeaks_plus))
            plt.xlabel('$B$ (T)', fontsize=20)            
            plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=20)  
            plt.xticks(fontsize=15)
            plt.yticks(fontsize=15)
            plt.gcf().subplots_adjust(bottom=0.15)
            plt.gcf().subplots_adjust(left=0.17)
            plt.title("Signál ze souboru "+str(counter_angle)+", úhel: "+'{:.1f}'.format(load_angle(counter_angle))+"°", fontsize=20)
        
        if plot_angle_sum == True:
            x_hpeaks_plus_sum = []
            x_lpeaks_minus_sum = []
            largest_val_hpeaks_sum = max(angle_hPeaks.shape)
            largest_val_lpeaks_sum = max(angle_lPeaks.shape)
            for y_hvars in range(largest_val_hpeaks_sum):  
                x_hpeaks_actual_sum = angle_hPeaks[:,0][y_hvars]
                if x_hpeaks_actual_sum >= 0:
                    ax6.plot(x_hpeaks_actual_sum, np.full(x_hpeaks_actual_sum.shape, load_angle(counter_angle)), 'ro')
                    x_hpeaks_plus_sum.append(x_hpeaks_actual_sum)
            for y_lvars in range(largest_val_lpeaks_sum):  
                x_lpeaks_actual_sum = angle_lPeaks[:,0][y_lvars]
                if x_lpeaks_actual_sum <= 0:
                    ax6.plot(x_lpeaks_actual_sum, np.full(x_lpeaks_actual_sum.shape, load_angle(counter_angle)), 'go')
                    x_lpeaks_minus_sum.append(x_lpeaks_actual_sum)
            x_lpeaks_minus_sum=np.sort(x_lpeaks_minus_sum)
            x_hpeaks_plus_sum=np.sort(x_hpeaks_plus_sum)
            
            if custom_peak_delete == True:
                if counter_angle >= 14:
                    x_hpeaks_plus_sum=np.delete(x_hpeaks_plus_sum, 2)  
                    
            ax2.plot(x_hpeaks_plus_sum, np.full(x_hpeaks_plus_sum.shape, load_angle(counter_angle)), 'ro')
            ax2.plot(x_lpeaks_minus_sum, np.full(x_lpeaks_minus_sum.shape, load_angle(counter_angle)), 'ro')
            
            cosine = np.cos((np.pi/180)*load_angle(counter_angle))
            if plot_Bperpendicular == True:
                 ax3.plot(x_hpeaks_plus_sum*cosine, np.full(x_hpeaks_plus_sum.shape, load_angle(counter_angle)), 'bo')
                 ax3.plot(x_lpeaks_minus_sum*cosine, np.full(x_lpeaks_minus_sum.shape, load_angle(counter_angle)), 'bo')
            print("Negative magnetic field peaks for angle: "+'{:.1f}'.format(load_angle(counter_angle))+"°: "+str(x_lpeaks_minus_sum))
            print("Positive magnetic field peaks for angle: "+'{:.1f}'.format(load_angle(counter_angle))+"°: "+str(x_hpeaks_plus_sum))
            
            datalist_angle.append(x_hpeaks_plus_sum)
            sizelist_angle.append(x_hpeaks_plus_sum.shape)
            list_y.append(np.full((last_measurement_data - 1), cosine))
            if plot_volt_dependance == True:
                ax4.plot((1/x_hpeaks_plus_sum), np.full(x_hpeaks_plus_sum.shape, cosine), 'ro')
                
        counter_angle += 1
        
    if (plot_angle_sum & plot_volt_dependance) == True:
        print("List of all detected peaks")
        print(datalist_angle)
        print("List of shapes of these peaks")
        print(sizelist_angle)
        print("Largest value on previous list")
        largest_val_angle=max(sizelist_angle)
        print(largest_val_angle[0])
        counter_angle = counter_angle_original
        apk1 = np.array([])
        apk2 = np.array([])
        apk3 = np.array([])  
        apk4 = np.array([])
        apk5 = np.array([])
        apk6 = np.array([])
        for sup_x in range(17):
            for sup_y in range(largest_val_angle[0]):            
                if sup_y == 0:
                    apk1 = np.append(apk1,datalist_angle[sup_x][sup_y])
                if sup_y == 1:
                    apk2 = np.append(apk2,datalist_angle[sup_x][sup_y])
                if sup_y == 2:
                    apk3 = np.append(apk3,datalist_angle[sup_x][sup_y]) 
                if sup_y == 3:
                    apk4 = np.append(apk4,datalist_angle[sup_x][sup_y]) 
                if sup_y == 4:
                    apk5 = np.append(apk5,datalist_angle[sup_x][sup_y]) 
                if sup_y == 5:
                    apk6 = np.append(apk6,datalist_angle[sup_x][sup_y]) 
                    
            pk_angle = np.append(pk_angle,load_angle(counter_angle))        
            counter_angle += 1
            
        ax4.plot(1/apk1, list_y, 'go')
        ax4.plot(1/apk2, list_y, 'ko')
        ax4.plot(1/apk3, list_y, 'bo')
        ax4.plot(1/apk4, list_y, 'yo')
        ax4.plot(1/apk5, list_y, 'co')
        ax4.plot(1/apk6, list_y, 'mo')
        beq1 =  np.polyfit((1/apk1), list_y, 1)
        beq2 =  np.polyfit((1/apk2), list_y, 1)
        beq3 =  np.polyfit((1/apk3), list_y, 1)
        beq4 =  np.polyfit((1/apk4), list_y, 1)
        beq5 =  np.polyfit((1/apk5), list_y, 1)
        beq6 =  np.polyfit((1/apk6), list_y, 1)
        ax4.plot(1/apk1, beq1[0]*(1/apk1)+beq1[1],'g')
        ax4.plot(1/apk2, beq2[0]*(1/apk2)+beq2[1],'k')
        ax4.plot(1/apk3, beq3[0]*(1/apk3)+beq3[1],'b')
        ax4.plot(1/apk4, beq4[0]*(1/apk4)+beq4[1],'y')
        ax4.plot(1/apk5, beq5[0]*(1/apk5)+beq5[1],'c')
        ax4.plot(1/apk6, beq6[0]*(1/apk6)+beq6[1],'m')
            
        print("Equation of first peak set (green): " +str(beq1[0][0])+" * x + "+str(beq1[1][0]))
        print("Equation of second peak set (black): " +str(beq2[0][0])+" * x + "+str(beq2[1][0]))
        print("Equation of third peak set (blue): " +str(beq3[0][0])+" * x + "+str(beq3[1][0]))
        print("Equation of fourth peak set (yellow): " +str(beq4[0][0])+" * x + "+str(beq4[1][0]))
        print("Equation of fifth peak set (cyan): " +str(beq5[0][0])+" * x + "+str(beq5[1][0]))
        print("Equation of sixth peak set (magenta): " +str(beq6[0][0])+" * x + "+str(beq6[1][0]))
       
if run_angle_comparison == True:
    fig5, ax5=plt.subplots()
    plt.axis([Fr_comp_min, Fr_comp_max, -8.50, 59.3])
    plt.xlabel('$B$ (T)', fontsize=25)            
    plt.ylabel('$Úhel$ (°)', fontsize=25)  
    plt.title("Úhel od "+'{:.1f}'.format(load_angle(counter_angle_original))+"° do "+'{:.1f}'.format(load_angle(last_measurement_data))+"°, vývoj s B", fontsize=25)
    plt.xticks(fontsize=25)
    plt.yticks(fontsize=25)
    
    fig6, ax6=plt.subplots()
    plt.axis([-0.275, 0.275, -8.50, 59.3])
    plt.xlabel(r'$B_{\bot}$ (T)', fontsize=25)            
    plt.ylabel('$Úhel$ (°)', fontsize=25)  
    plt.title("Úhel od "+'{:.1f}'.format(load_angle(counter_angle_original))+"° do "+'{:.1f}'.format(load_angle(last_measurement_data))+"°, vývoj s kolmou složkou B", fontsize=25)
    plt.xticks(fontsize=25)
    plt.yticks(fontsize=25)
        
    counter_angle = counter_angle_original
    while counter_angle <= last_measurement_data:
        loaded_angle_data = load_angle_data(counter_angle)
        AngleData = np.loadtxt(open(loaded_angle_data,'r'),skiprows=30) 
        loaded_offset = load_offset(counter_angle)
        
        B_ang = (AngleData[:,0])+loaded_offset
        ds_ang = AngleData[:,1]
        B_ang_det = B_ang
        ds_ang_det = ds_ang
        B_ang = B_ang[np.logical_and(-Fr_val<B_ang,B_ang<Fr_val)]
        B_ang_det = B_ang_det[np.logical_and(Fr_comp_min<B_ang_det,B_ang_det<Fr_comp_max)]
        ds_ang = ds_ang[np.logical_and(-Fr_val<B_ang,B_ang<Fr_val)]
        ds_ang_det = ds_ang[np.logical_and(Fr_comp_min<B_ang,B_ang<Fr_comp_max)]
        ds_smooth = ds_ang_det
        
        '''Measurement data normalization'''
        angle_normalized = []
        angle_max = abs(max(ds_smooth, key=abs))
        for pvar in ds_smooth:
            angle_normalized.append((((pvar/angle_max)*multiplicator_angle_comparison)+load_angle(counter_angle)))
         
        cosine = np.cos((np.pi/180)*load_angle(counter_angle))
        if (plot_angle_sum & plot_volt_dependance) == True:
            apk1_val = (np.nonzero(B_ang_det[:,None] == apk1))[0]
            apk2_val = (np.nonzero(B_ang_det[:,None] == apk2))[0]
            apk3_val = (np.nonzero(B_ang_det[:,None] == apk3))[0]
            apk4_val = (np.nonzero(B_ang_det[:,None] == apk4))[0]
            apk5_val = (np.nonzero(B_ang_det[:,None] == apk5))[0]
            apk6_val = (np.nonzero(B_ang_det[:,None] == apk6))[0]
            if counter_angle == 3:
                apk1_val = apk1_val[0]
                apk2_val = apk2_val[0]
                apk3_val = apk3_val[0]
                apk4_val = apk4_val[1]
                apk5_val = apk5_val[0]
                apk6_val = apk6_val[0]
            elif counter_angle == 4:
                apk1_val = apk1_val[1]
                apk2_val = apk2_val[0]
                apk3_val = apk3_val[0]
                apk4_val = apk4_val[0]
                apk5_val = apk5_val[0]
                apk6_val = apk6_val[0]
            elif counter_angle == 8:
                apk1_val = apk1_val[0]
                apk2_val = apk2_val[0]
                apk3_val = apk3_val[1]
                apk4_val = apk4_val[0]
                apk5_val = apk5_val[1]
                apk6_val = apk6_val[0]
            elif counter_angle == 9:
                apk1_val = apk1_val[0]
                apk2_val = apk2_val[0]
                apk3_val = apk3_val[0]
                apk4_val = apk4_val[0]
                apk5_val = apk5_val[1]
                apk6_val = apk6_val[0]
            else:
                apk1_val = apk1_val[0]
                apk2_val = apk2_val[0]
                apk3_val = apk3_val[0]
                apk4_val = apk4_val[0]
                apk5_val = apk5_val[0]
                apk6_val = apk6_val[0]
                
            ax6.plot(B_ang_det[apk1_val]*cosine, angle_normalized[apk1_val], 'kx')
            ax6.plot(B_ang_det[apk2_val]*cosine, angle_normalized[apk2_val], 'kx')
            ax6.plot(B_ang_det[apk3_val]*cosine, angle_normalized[apk3_val], 'kx')
            ax6.plot(B_ang_det[apk4_val]*cosine, angle_normalized[apk4_val], 'kx')
            ax6.plot(B_ang_det[apk5_val]*cosine, angle_normalized[apk5_val], 'kx')
            ax6.plot(B_ang_det[apk6_val]*cosine, angle_normalized[apk6_val], 'kx')
            
        ax5.plot(B_ang_det, angle_normalized, 'b')
        ax6.plot(B_ang_det*cosine, angle_normalized, 'r')
        
        if run_overlap == True:
            Ef = Ef_start
            file_title = 'Ef_Sim/Ef_' +  '{:.1f}'.format(Ef*1000) + 'meV.npz'
            Loaded_data = np.load(file_title)
                
            '''Simulation data normalization'''
            ds_xx_normalized = []
            dsval_max = abs(max(Loaded_data['ds_xx'], key=abs))
            for dsvar in Loaded_data['ds_xx']:
                ds_xx_normalized.append(((dsvar/dsval_max)*2.2)+load_angle(counter_angle))
                
            minus_B = -1*Loaded_data['B'][::-1]
            minus_ds = ds_xx_normalized[::-1]
            plus_B = Loaded_data['B']
            plus_ds = ds_xx_normalized
            ax5.plot(minus_B/cosine, minus_ds, 'k', linewidth=0.75)
            ax5.plot(plus_B/cosine, plus_ds, 'k', linewidth=0.75)

            del Loaded_data.f
            Loaded_data.close()
        counter_angle += 1
     
plt.show()
