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

def load_data(dummy_counter): #Containt measurement data locations
    switcher={
        1: "2008-09-29 Natural Graphite - SdH temps/01_100K-285GHz.dat",
        2: "2008-09-29 Natural Graphite - SdH temps/46_15K-285GHz.dat",
        3: "2008-09-29 Natural Graphite - SdH temps/35_07K-285GHz.dat",
        4: "2008-09-29 Natural Graphite - SdH temps/45_05K-285GHz.dat"      
        }
    return switcher.get(dummy_counter, "Error: File does not exist")
f = 285 #Measurement frequency
T = [100, 15, 7.5, 5] #Temperatures during measurement
det_B_min = 1 #Minimal magnetic field for plotting
det_B_max = 7 #Maximal magnetic field for plotting
color = ['r', 'b', 'g', 'm']
counter = 2 #From which file to start
smooth_const = 15 #constant for data smoothing

datalist_s = []
fig,ax1=plt.subplots()
plt.xlabel(r'$B$ (T)', fontsize=25)
plt.ylabel(r'$d\sigma_{xx}/dB$ (arb.u.)', fontsize=25)
plt.title(r'$f$='+str(f)+' GHz, srovnání dat derivace optické konduktivity', fontsize=25)
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)

fig,ax2=plt.subplots()
plt.xlabel(r'$B$ (T)', fontsize=25)
plt.ylabel(r'$\sigma_{xx}$ (arb.u.)', fontsize=25)
plt.title(r'$f$ = '+str(f)+' GHz, srovnání dat optické konduktivity', fontsize=25)
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)

while counter <= 4:
        '''Data load'''
        loaded_data = load_data(counter)
        DataFile = np.loadtxt(open(loaded_data,'r'),skiprows=30) 
        B_exp=DataFile[:,0]
        ds_exp = DataFile[:,1]
        if counter == 4:
            B_exp = B_exp[::-1]
            ds_exp = ds_exp[::-1]
        B_exp_det = B_exp
        ds_exp_det = ds_exp
        B_exp = B_exp[np.logical_and(-16<B_exp,B_exp<16)]
        B_exp_det = B_exp_det[np.logical_and(det_B_min<B_exp_det,B_exp_det<det_B_max)]
        ds_exp = ds_exp[np.logical_and(-16<B_exp,B_exp<16)]
        ds_exp_det = ds_exp[np.logical_and(det_B_min<B_exp,B_exp<det_B_max)]
        s_exp_det = integrate.cumtrapz(ds_exp_det,B_exp_det)
        B_exp_det2 = B_exp_det[1:]
        ds_smooth = ff.smooth(ds_exp_det, smooth_const)
        s_smooth = s_exp_det
        Integration_normalized = s_smooth
        Derivation_normalized = ds_smooth
        '''Baseline subtraction'''
        lin_params=np.polyfit(B_exp_det2,Integration_normalized,1)
        base=ff.smooth(Integration_normalized, 501)
        Integration_normalized = Integration_normalized - base
        
        '''Peak Detect'''
        peaks = pk.peakdetect(Integration_normalized, B_exp_det2, lookahead=10, delta=1e-9) 
        lPeaks = np.array(peaks[1])
        datalist_s.append(lPeaks[:,1][0])
        hPeaks = np.array(peaks[0])
        datalist_s.append(hPeaks[:,1][0])
        peaks_all = pk.peakdetect(Integration_normalized, B_exp_det2, lookahead=5, delta=5e-7) 
        hPeaks_all = np.array(peaks_all[0])      
        ax1.plot(B_exp_det,Derivation_normalized, color[counter-1], label = (str(T[counter-1])+' K'))        
        ax2.plot(B_exp_det2,Integration_normalized, color[counter-1], label = (str(T[counter-1])+' K'))
        
        fig,ax=plt.subplots()
        ax.plot(B_exp_det2,Integration_normalized+base, label=r'$\sigma_{xx}$')
        ax.plot(B_exp_det2,base, label='Baseline')
        ax.plot(B_exp_det2, Integration_normalized, label = r'$\sigma_{xx}$ background-corrected')
        #ax2.plot(hPeaks_all[:,0], hPeaks_all[:,1], 'ro')
        plt.xlabel('$B$ (T)', fontsize=16)
        plt.ylabel('Microwave absorption (arb. u.)', fontsize=16)
        plt.title(r'$f$='+str(f)+' GHz, $T$='+str(T[counter-1])+' K', fontsize=16)
        plt.xticks(fontsize=13)
        plt.yticks(fontsize=13)
        plt.legend(fontsize=14)
        
        fig,ax5=plt.subplots()
        counter_aux = 0
        sizelist = []
        while counter_aux < len(hPeaks_all):
            sizelist.append(counter_aux+1)
            counter_aux += 1
        params=np.polyfit(1/hPeaks_all[:,0],sizelist,1)
        line = params[0]*(1/hPeaks_all[:,0]) + params[1]
        ax5.plot(1/hPeaks_all[:,0], line, color[counter-1])
        ax5.plot(1/hPeaks_all[:,0], sizelist, 'ro')
        plt.xlabel('$1/B$ (T)', fontsize=16)
        plt.ylabel('Peak number', fontsize=16)
        plt.title(r'Peaks position, $T$='+str(T[counter-1])+' K', fontsize=16)
        print("Function for T="+str(T[counter-1])+"K: "+str(params[0])+"*x+"+str(params[1]))
        counter += 1
   
print("Peaks list for temps:")
print(datalist_s)
min_temp = datalist_s[2]
print(min_temp)
ax1.legend(fontsize=20)
ampl=np.array([8.4E-5,6.18E-5,2.41E-5])
Temp=np.array([5,7.5,15])
Temp2 = np.arange(3,30,0.1)

def Dingle(T,T0,Amstar): #Fitting function calculation
    return T/T0 * np.sinh(Amstar)/np.sinh(Amstar*T/T0)
Tparams,desv_est = optimize.curve_fit(lambda T,Amstar: Dingle(T,Temp[0] ,Amstar), Temp,ampl/ampl[0])
    
fig,ax1=plt.subplots()
ax1.plot(Temp,ampl/ampl[0], 'ro')
ax1.plot(Temp2,  Dingle(Temp2, Temp[0],Tparams[0]))
plt.xlabel(r'$T$ (K)', fontsize=16)
plt.ylabel(r'$(d\sigma_{xx}/dB)_{T}/(d\sigma_{xx}/dB)_{Tmin}$ (-)', fontsize=16)
plt.title('Velikost amplitudy v závislosti na teplotě', fontsize=16)
A = (2*(pi**2)*kb*5)/((h/2*pi)*e*2)
print(Tparams)
print(desv_est)
print(A)
result=(Tparams[0]/A)/m_e
print(result)
fig,ax3=plt.subplots()
ax3.plot(T[1], datalist_s[0]/min_temp, "bo", label = '15 K')
ax3.plot(T[2], datalist_s[1]/min_temp, "go", label = '7,5 K')
ax3.plot(T[3], datalist_s[2]/min_temp, "mo", label = '5 K')
plt.xlabel(r'$T$ (K)', fontsize=16)
plt.ylabel(r'T/Tmin (arb.u.)', fontsize=16)
plt.title('SdH amplitude height in dependance on temperature', fontsize=16)
ax3.legend(fontsize=14)
ax2.legend(fontsize=20)

'''FOURIER TRANSFORM CALCULATION'''
fileName="Experimental data/004_Graphite_Suhner.dat"
Data = np.loadtxt(open(fileName,'r'), skiprows=30)
B_exp=Data[:,0]
ds_exp=Data[:,1] #Derivative of conductivity
dsy_exp=Data[:,2]
s_exp = integrate.cumtrapz(ds_exp,B_exp) #conductivity, integrating reduces the length of the array by one
B_exp2= B_exp[1:] #reducing the magnetic field by one entry

'''Selecting the magnetic field range of the date'''
arr1= np.flatnonzero(np.logical_and(B_exp>0.68,B_exp<1))
arr2=arr1
B_p=B_exp[arr1]
ds_p=ds_exp[arr1]
B_p2=B_exp2[arr2]
s_p=s_exp[arr2]*1E6 #1E6 to get the signal in the order of 1E1
ds_smooth=ds_p
s_smooth=s_p

'''subtracting the baseline '''
#base by linear fit
lin_params=np.polyfit(B_p2,s_smooth,1)
base= lin_params[1] + lin_params[0]*B_p2
s_smooth = s_smooth - base

''' Getting the local maxima '''
B_max,i_max=ff.Bminima(-1*s_smooth,B_p2,Bwindow=0.01) #multiplying by -1 to get the maxima.
print('Bmax ', B_max)

'''Calculating the Fourier transform'''
#Algorithm made by Jorge. Not good results as there is no windowing of the data
X,Y,x2,y2=ff.fourier_SdH(s_smooth,B_p2) #x2 = 1/B. y2= rearanged s_smooth

###
#calculating the Fourier transf. with other algorithm
# found in page https://dsp.stackexchange.com/questions/24680/why-is-the-magnitude-of-the-second-component-of-my-fft-spectrum-always-the-large?rq=1 
#for zero padding the signal
#exponent for the next power of two higher than the length of the signal. for instace, len(y2)==10, pad=4, 2**(pad)=16. padded y2:  [y2[0],...,y2[9],0,0,0,0,0,0,0]
pad= np.log10(len(y2))/np.log10(2) 
pad=int(pad)+1 
# using Hann window and padding
freq, spectrum = signal.periodogram(y2, window="hann", fs=1/(x2[1]-x2[0]), nfft=2**pad) 

#Physical parameters for the Fourier transform
mc0=0.06*m_e #effective mass of the carriers at the H point of graphite, in kg
f=283.2E9 #Frequency of the measurement in Hz
w=2*np.pi*f # Microwave frequency in Hz
wc0= e*B_p/mc0 #cyclotron frequency w.r.t the magnetic field. in Hz

fig,ax=plt.subplots()
ax.plot(B_exp,ds_exp*3333)
plt.xlabel('$B$ (T)', fontsize=20)
plt.ylabel(r'$d\sigma_{xx}/dB$ (arb.u.)', fontsize=20)
plt.xticks(fontsize=13)
plt.yticks(np.arange(-1, 1.1, step=0.2), fontsize=13)
plt.axis([-1.1,1.1, -1.1, 1.1])
plt.title('$f=283.2$ GHz, $T=5$ K', fontsize=20)
plt.gcf().subplots_adjust(left=0.15)
plt.gcf().subplots_adjust(bottom=0.14)

fig,ax=plt.subplots()
ax.plot(B_p,ds_smooth*3333)
plt.xlabel('$B$ (T)', fontsize=20)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.axis([0.66,1.02, -0.4, 0.4])
plt.ylabel('$d\sigma_{xx}/dB$ (arb.u.)', fontsize=20)
plt.title("Vybraná oblast dat SdHO",fontsize=20)
plt.gcf().subplots_adjust(left=0.15)
plt.gcf().subplots_adjust(bottom=0.14)

fig,ax=plt.subplots()
ax.plot(B_p2,s_smooth+base, label=r'$\sigma_{xx}$')
ax.plot(B_p2,base, label='Základní linie')
ax.plot(B_p2,s_smooth, label=r'$\sigma_{xx}$ s korekcí')
plt.xlabel('$B$ (T)', fontsize=20)
plt.ylabel('$\sigma_{xx}$ (arb.u.)', fontsize=20)
plt.title("Numerická integrace vybrané oblasti",fontsize=20)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.legend(fontsize=15)
plt.gcf().subplots_adjust(left=0.15)
plt.gcf().subplots_adjust(bottom=0.14)

fig,ax=plt.subplots()
ax.plot(x2,y2)
plt.xlabel(r'$1/B$ (T$^{-1}$)', fontsize=20)
plt.ylabel('$\sigma_{xx}$ (arb.u.)', fontsize=20)
plt.title("Převrácená hodnota $B$  pro $\sigma_{xx}$",fontsize=20)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.axis([0.95,1.5, -0.5, 0.5])
plt.gcf().subplots_adjust(left=0.15)
plt.gcf().subplots_adjust(bottom=0.14)

#FFT using PSD
fig,ax=plt.subplots()
ax.plot(freq,np.abs(spectrum)*645)
plt.axis([0,40,0,750*np.max(np.abs(spectrum))])
peak_x, peak_y= freq[np.argmax(np.abs(spectrum))], 645*np.max(np.abs(spectrum))
dx= ax.get_xlim()[1]-ax.get_xlim()[0] #length of the plot in the x axis
ann_text= r'$B_F = $' + '{:.2f}'.format(peak_x) + ' T'
plt.annotate(ann_text, xy=(peak_x,peak_y), xytext=(peak_x+dx*0.1 ,peak_y),arrowprops=dict(arrowstyle="-|>",facecolor='black',relpos=(0, 0.5)), fontsize=13)
plt.title('Fourierova Transformace', fontsize=15)
plt.xlabel('Frekvence FT (T)', fontsize=15)
plt.ylabel('Amplituda FT (arb. u.)', fontsize=15)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)

plt.show()
